#include "mex.h"
#include <math.h>

static void scalarxarray(double sca, double arr[], int N)
{
    int i;
    for(i=0;i<N;i++)
        arr[i]=sca*arr[i];

}

static void neo(double* asset,int Npath,int Dtime)
{
    int i,j;
    for (j=0;j<Dtime;j++)
            for(i=0;i<Npath;i++)
                asset[j*Npath+i]=0.0;
}

static void pricex1(double* asset, double* eta, double* etasquare, double *etaeta, double *X1,int Npath, int betaK)
{
    int i,j,k,iter, counter=0,column=betaK*(betaK-1)/2;
    double currentx;
    //double* asset=malloc(Npath*sizeof(double));
    
    
    
//     for iter=1:51
// asset=zeros(Npath,1);
// for ii=2:betaK+1
//     asset=asset+X(iter,ii)*eta(:,ii-1);
// end
// for ii=betaK+2:2*betaK+1
//     asset=asset+X(iter,ii)*etasquare(:,ii-betaK-1);
// end
// for ii=2*betaK+2:2*betaK+1+betaK*(betaK-1)/2
//     asset=asset+X(iter,ii)*etaeta(:,ii-2*betaK-1);
// end
// asset=asset+X(iter,1);
// 
// price(iter)=sum(abs(asset))/Npath;
// end
    
 
        for (j=1;j<betaK+1;j++)
        {
            for(i=0;i<Npath;i++)
                asset[i]=asset[i]+X1[j]*eta[(j-1)*Npath+i];
        }
        for (j=betaK+1;j<2*betaK+1;j++)
        {
            for(i=0;i<Npath;i++)
                asset[i]=asset[i]+X1[j]*etasquare[(j-betaK-1)*Npath+i];
        }
        for (j=2*betaK+1;j<2*betaK+1+betaK*(betaK-1)/2;j++)
        {
            for (i=0;i<Npath;i++)
                asset[i]=asset[i]+X1[j]*etaeta[(j-2*betaK-1)*Npath+i];
        }
        for (i=0;i<Npath;i++)
            asset[i]=asset[i]+X1[0];
}


void mexFunction(int nlhs,mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
        double *asset, *eta,*etasquare,*etaeta, *X1;
        int Npath,betaK,Dtime;
        int i;
        
        Npath=mxGetM(prhs[0]); 
        betaK=mxGetN(prhs[0]);
//        Dtime=mxGetM(prhs[3]);
        
        plhs[0]=mxCreateDoubleMatrix(Npath,1,mxREAL);
                
//        plhs[1]=mxCreateDoubleMatrix(Dtime,1,mxREAL);
        
        asset=mxGetPr(plhs[0]);
        
        eta=mxGetPr(prhs[0]);
        etasquare=mxGetPr(prhs[1]);
        etaeta=mxGetPr(prhs[2]);
        
        X1=mxGetPr(prhs[3]);
        
        
        
        //neo(asset,Npath,Dtime);
        pricex1(asset,eta,etasquare,etaeta,X1,Npath,betaK);


}
